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EXCUTIVE SUMMARY 


A fully 3-D numerical model is developed to represent magnetic damping of complex 
fluid flow, heat transfer and electromagnetic field distributions in a melt cavity. The model is 
developed based on our in-house finite element code for the fluid flow, heat transfer and 
electromagnetic field calculations. The computer code has been tested against benchmark test 
problems that are solved by other commercial codes as well as analytical solutions whenever 
available. The numerical model is tested against numerical and experimental results for water 
reported in literature. With the model so tested, various numerical simulations are carried out for 
the Sn-35.5% Pb melt convection and temperature distribution in a cylindrical cavity with and 
without the presence of a transverse magnetic field. Numerical results show that magnetic 
damping can be effectively applied to reduce turbulence and flow levels in the melt undergoing 
solidification and over a certain threshold value a higher magnetic field resulted in a higher 
velocity reduction. It is found also that for a fully 3-D representation of the magnetic damping 
effects, the electric field induced in the melt by the applied DC magnetic field does not vanish, as 
some researchers suggested, and must be included even for molten metal and semiconductors. 
Also, for the study of the melt flow instability, a long enough time has to be applied to ensure the 
final fluid flow recirculation pattern. Moreover, our numerical results suggested that there seems 
to exist a threshold value of applied magnetic field, above which magnetic damping becomes 
possible and below which the convection in the melt is actually enhanced. Because of the limited 
financial resource allocated for the project, we are unable to carry out extensive study on this 
effect, which should warrant further theoretical and experimental study. In that endeavor, the 
developed numerical model should be very useful; and the model should serve as a useful tool 
for exploring necessary design parameters for planning magnetic damping experiments and 
interpreting the experimental results. 
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1. INTRODUCTION 

Recently, NASA Lewis has conducted a set of Sn-Pb alloy solidification experiments in 
High Temperature Directional Solidification Furnace under the influence of a transverse 
magnetic field to investigate magnetic damping effects (Song, Tewari and de Groh, 1996). The 
schematical representation of their system is sketched below in Figure 1. 
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Figure 1. Schematic representation of the solidification processes at NASA Lewis. 


With the experimental setup, they studied fluid flow, thermal transport and 
macrosegregation during solidification of Pb-Sn alloys with and without an externally imposed 
horizontal magnetic field. Temperature fluctuations were observed in the liquid pool during the 
solidification processes. It was found that these thermal fluctuations evolve from none to cyclic 
to time periodic having multiharmonics and finally to random as the Rayleigh number of the 
system increases. The application of the magnetic field decreases convection level in the liquid 
pool, but even high magnetic field (0.45 T) was unable to eliminate the convection. Application 
of the magnetic fields also resulted in changing temperature fluctuations from sinusoidal to 
pulsed wave behavior. These experimental measurements are very useful for providing first 
hand evidence on the effectiveness of magnetic damping. However, accompanying theoretical 
analyses are needed to determine quantitatively how the thermal field is affected by a 
magnetically damped fluid flow field during solidification. 

The objective of this small grant is to initiate the development of first principle 
mathematical models for magnetically damped fluid flow and heat transfer in solidifying liquids. 
Our intention is two fold: (1) to acquire a basic understanding of magnetically damped flow 
phenomena in the aforementioned system by studying a simplified system such as a 3-D 
cylindrical cavity without considering the complex phase change phenomena associated with 
solidification; and (2) to establish a capability for 3-D simulation of magnetically damped flows 
with which further studies on magnetic damping effects on melt flow and solidification processes 
can be conducted. Following this direction, we have developed a full 3-D finite element model, 
based on our in-house fluid flow and heat transfer and electromagnetic model capabilities, for 
magnetic damping of fluid flow and thermal phenomena in a 3-D cylindrical cavity. The 
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diameter of the cavity is the same as used in the Lewis experiments (Song et al., 1996), but the 
length is varied to simulate the different stages of solidification. Sn-Pb melt is used and the 
properties of the melt are provided by H. C. de Groh at NASA Lewis, along with other operating 
parameters. Both pseudo steady state and transient numerical models were developed. Some 
simulations were conducted with the model developed to both demonstrate the capability of the 
model and also to gain some basic understanding of magnetic damping effect during melt 
processing of relevance to NASA Lewis experiments. The modeling strategy, numerical results 
and findings and the future directions on the work are described below. 

2. THE 3-D FINITE ELEMENT MODEL 

The mathematical formulation of the magnetic damped fluid flow and heat transfer 
problem is given in Appendix I. Because of the limited funds available, we focus our attention 
on the model capability development and basic understanding, rather than a comprehensive 
study, of the magnetic damping effects on fluid flow and heat transfer phenomena in the Sn-Pb 
melts during solidification. Thus, we consider natural convection and magnetic damping 
phenomena in a 3-D cylindrical cavity with its side thermally insulated and top and bottom 
surfaces fixed at lower and high temperatures. Solidification and other complicating factors such 
as crystal growth rate, etc., are not considered in the model to simplify the calculations. 

Originally the models were planned to be developed using FIDAP, a commercial finite 
element code for flow and thermal computations, and simulations were to be carried out using 
the Supercomputers at the University of Illinois. However, two difficulties encountered during 
the course of the model development with FEDAP. First several testing runs with FIDAP for 
flows without applied magnetic fields were made on the Supercomputer and it was found that it 
took too long to obtain the results. The major time was spent on waiting and queue list, even 
though the highest priority was employed. For example, for a 3-D mesh of 12 x 5 x 10 8-node 
elements used for the simulation of natural convection in a cylinder of our interest, one transient 
run that requires about 200 time steps took about 4 CPU hours; but the real turnout time was 
about 4 to 5 days, because of the heavy usage of the machine. Second, some time during the 
model development when we carefully studied FIDAP’s formulation of magnetic damping 
problems as stated in their manual (p. 2-32, FIDAP Theoretical Manual), we found that the 
current density resulting from the interaction of the magnetic field with the fluid flow field seems 
to be calculated incorrectly. Since the source code is not available, we were not able to confirm 
whether this is a typographical error or the code calculates the problem as stated (usually it does 
and should also). Because of the limited computing sites where FIDAP is available and because 
there is potentially an error in their problem formulation, we started to develop our own 3-D 
finite element code for the planned numerical modeling effort. Fortunately, at that time, our 2-D 
finite element model was completed and was in use for g-jitter induced fluid flow, heat and mass 
transfer calculations. Our 3-D code development was then just an extension of the 2-D finite 
element code. The finite element formulation of the magnetic damping problems based on 
which our code was developed is provided in Appendix II. The development of our own code 
entailed program coding, debugging and extensive testing of the code against analytical solutions 
and calculations from other codes such as FEDAP and FLOW3D for benchmark type problems. 
Some of testing has been discussed in one of our recent publications (Pan, Li and deGroh, 1997). 
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It should be noted that having our own code gives us tremendous flexibility, although the 
development process itself is a bit painful in that various cases, including those not of direct 
interest to our model development, must be examined to ensure the correctness of the code. 
With our own code, we now are able to load the program onto any platforms we can get an 
access to and thus many computational cases can be run simultaneously to increase the 
computational efficiency. Many calculations were carried out on the LSU EBM-6000 RISC 
cluster machines, where several runs can be made simultaneously. 

In Fall 1997, the PI left LSU to join the School of Mechanical and Materials Engineering 
at Washington State University. The model was imported from Louisiana State University to 
Washington State University and additional cases were run on an Onix machine, which is 
capable of supercomputing. For the results presented below, a typical run of 2000 time steps 
leading to a steady state would take 20 CPU hours for a mesh size of 17x6x10 8-node elements 
on the Onix machine. This change in environment required re-setting up the local computing 
environment for the Pi’s laboratory, including purchase of PCs, printers and tuning software, 
etc., which are required for printing out the calculated results and preparing this report. 

3. NUMERICAL RESULTS AND FINDINGS 

The developed numerical model is used to carry out simulations for magnetically damped 
flows in Sn-Pb melts for various operating conditions. The parameters for the calculations are 
given in Tables 1 and 2. To test the mesh dependency of our calculations, numerical 
experiments were carried out using different mesh sizes or different number of elements. 
Basically, different meshes were used until the results no longer show the mesh dependency. 
The mesh, as shown in Figure A 2, is the one of the meshes used for the calculations for cavities 
with small h/d ratio. More meshes are placed in the vertical direction when a larger h/d ratio 
cavity is considered. Unless stated otherwise, all the calculations for magnetic damping assumed 
that the magnetic field is applied in the x-direction (see Figure Al). For the results presented in 
the figures below, the time and velocity are normalized by the time factor (t) and the velocity 
factor (Uo) given in Table 2 and Rayliegh numbers and other conditions are based on the stability 
diagram given by Muller (1988). 

(a) Test against published results for known conditions 

Before reporting the results on Sn-35.5% Pb melt flows, we would like to discuss some 
of testing cases with water for which published work is available. Muller’s research group had 
published results on flow instability of water in a cylinder heated from below for various aspect 
ratios and also conducted experiments to very their calculations (Muller et al., 1984). Figures 2 
and 3 depict the comparison between our computed results and the results from a finite 
difference model and corresponding water model developed by Muller’s group. The parameters 
used for the calculations are given in Table 3. It is clear that for both symmetric and un- 
symmetric flow patterns the results computed from our 3-D finite element model compared very 
well with those reported by Muller’s group, indicating that our numerical model can provide 
reasonably good prediction of fluid flow and heat transfer behavior of the fluids in the geometry 
under consideration. 
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(b) Melt convection in the cylindrical cavity without imposed magnetic Held 

Figures 4-6 show the evolution of the temperature and fluid flow fields for the Sn-35.5% 
Pb melt in the cylindrical cavity of h/d=0.5. Both 3-D and 2-D views (plane views) are 
presented. Note that the contour levels in Figure 4 correspond to those marked in Figure 6. It is 
apparent that both the fluid flow and temperature fields evolve starting from some arbitrary 
temperature pulse and eventually arrive at a steady state after 2000 time steps (corresponding to 
the real time of 3.05 seconds and dimensionless t=2.0). At the final steady state, a symmetric 
pattern is obtained with the flow moving downward at the center and upward from the side wall. 
It should be stressed that for such an unstable case a final stage of the flow pattern is strongly 
affected by the initial temperature pulse to induce the flow. When an appropriate temperature 
pulse is used, a flow pattern opposite to what is shown is also possible. This has been reported 
for water case and confirmed by our study as well (see Figure 2 above). 

The transient melt flow and temperature distribution for the Sn-35.5% Pb melt in a cavity 
of h/d=l is shown in Figures 7-9. In this case, a non-symmetric flow pattern is developed in the 
cylinder (see Figures 7(d), 8(d) and 9(g) and (h), which is caused by the non-symmetric pulse 
(see Figure 7(a)). Again, a long simulation time (time step=2000) is required to reach the steady 
state, afterwards the fields values are no longer changing. 

Similar transient calculations were also carried out for the melt column with a higher 
aspect ratio, h/d=3. The final steady state results are given in Figure 10. Clearly, the strong flow 
is developed near the bottom which is fixed at a higher temperature. The flow pattern is non- 
symmetric and the flow has a strong effect on the temperature field there as well. In the upper 
portion of the melt pool, the fluid is basically quiescent. 

(c) Magnetic damping effects on the melt flow in the cavity 

Figures 11-13 shows the transient development of the fluid flow and temperature 
distributions in the Sn-35.5% Pb melt in the cavity of h/d=0.5 with an applied magnetic field 
(B=0.05 T), with other conditions kept the same as those in Figure 4-6. Comparison of Figs 11- 
13 with Figs 4-6 clearly indicates the magnetic damping effects. In particular, when the system 
reaches the steady state, the fluid flow level is damped very substantially and temperature profile 
is basically flat. Note also that the steady state is reached within about 800 time steps. 

The magnetic damping effect on the melt flow in the Sn-35.5% Pb melt with an aspect 
ratio of h/d=l is shown in Figures 14 and 15. Apparently, with an increase in the applied 
magnetic field strength, magnetic damping becomes more effective, as expected. Figure 16 
shows the numerical results obtained for magnetic damping of melt flow in the Sn-35.5% Pb 
melt with an aspect ratio of h/d=3 in a transverse magnetic field (B=0.05T), where it can be seen 
that the applied magnetic field results in the reduction of the melt flow in the liquid pool. 
However, in comparison with Figure 10, the temperature field does not seem to change greatly. 

The higher magnetic field effect on the fluid flow and temperature distributions in the Sn- 
35.5% Pb melt for 3 different aspect ratios (h/d=0.5, h/d=l and h/d=3) is illustrated in Figures 
17-19, where B=0.5 T. Clearly, for all the cases, the damping effects are sufficiently strong that 
the temperature distribution in the system is basically controlled by diffusion, as is evident by the 
parallel temperature profiles in these figures. 
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One of the important reasons for applying the magnetic damping is to reduce the possible 
melt flow turbulence during solidification. The current model should help to provide some 
insight into the effectiveness of magnetic damping effects on the melt flow turbulence in the 
system. Figure 20 compares the effect of the different magnetic field strengths as applied to 
damp the turbulence in the system (with h/d=l, Ra=8 x 10 5 ). Inspection of the plots in the figure 
indicates that in the absence of the magnetic field (B=0) and with a small magnetic field present 
(B=0.005 T), the flow remains turbulent. With a higher magnetic field (B=0.05T), the 
turbulence is substantially reduced and melt flows in the laminar region. With an even higher 
magnetic field (B=0.5T), the turbulence is further reduced (the non-dimensionalized Umax is 
reduced to 1.57 from 172.5). Note that the velocity is plotted with the same magnitude scale for 
all the plots in Figure 20; the temperature contour clearly illustrates the effects of the melt flow 
still present in the system. 

Besides the results reported here, additional runs were made to investigate the effect of 
the direction of the applied magnetic field on the fluid flow and heat transfer effects. Results 
show that the magnetic field, when applied in the transverse direction provides the most effective 
damping effects. For symmetric flow pattern, the x- and y-direction magnetic field seems to 
produce similar results. For a non-symmetric flow pattern, the damping effect is stronger when 
the magnetic field is applied at the plane in which flow circulation is strongest. 

(d) Research findings that are different than those reported 

Our numerical simulations unvealed two important points that are either unclear in the 
literature or reported differently. For the Sn-35.5% Pb case of h/d=0.5, we found that 
calculations have to proceed for longer time than 600 (dimensionless) used by Muller’s group in 
order to resolve the ordered pattern. Figure 21 shows the evolution of the fluid flow and 
temperature profiles for the Sn-35.5% Pb melt at various time steps. According to Muller’s 
numerical study of molten melts with similar properties in which numerical simulations were 
carried out up to 600 time steps (undimensionalized the same as this study) (Muller, 1988), no 
symmetric flow pattern would exist for this condition. Our results compare well with theirs up to 
the time they terminated their calculations (time step=600) (Figure 21(a) shows the flow field at 
time step=500). However, the flow is actually still evolving and a longer run time is needed for 
this case to reach a steady state, as clearly indicated by the plotted results. For this particular 
case under consideration, the steady state is attained at about time step=2000 and by then a 
symmetric pattern is attained. It is noted here that Muller’s work might have used a different 
starting pulse because his flow pattern is just opposite to that reported here; nonetheless a steady 
state of symmetric pattern should prevail if a long run time is used. 

Another important finding resulted from this study, but unclear or misleading in the 
literature, is that for 3-D magnetic damping modeling studies, the induced electric field must be 
included in the current density and hence the Lorentz force calculations. Muller’s group recently 
did some 3-D simulations for semiconductor melt but without considering the electric field effect 
(Baumgartl et al., 1990). They argued, based on the cited references, that the conductivity of the 
molten metal and semiconductor is high and therefore the electric field is zero everywhere. Our 
study indicated that these references dealt with either 2-D or axi-symmetric case with a special 
arrangement of the magnetic field and for these special cases one can approve that the induced 
electric field is indeed zero everywhere in the melt. This approval, however, can not be 
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generalized to a fully 3-dimensional case. Indeed, our numerical calculations showed that for a 
fully 3-D model the induced electric field is not zero everywhere and hence must be included in 
the damping force calculations. Some of these calculated results are given in Figure 22. 

(e) Some unresolved issues that require further studies 

From our limited numerical experiments, we also found that a pure application of 
magnetic field does not necessarily decrease the fluid flow, as one would normally expect from 
simple Hartmann flow studies. On contrary, a DC magnetic filed may actually enhance the 
liquid convection in melts when a weak magnetic field is applied. A set of these results is given 
in Figure 23. There it is clear that with B=0.005T, the flow velocity is actually increased and 
damping effects are not obvious until the applied magnetic field is about 0.016 T. This 
phenomenon seems to occur in all other cases we tested. Our studies indicate that this flow 
augumentation seems to be dependent on the combination of the Rayleigh number of the system 
and the strength of the applied magnetic field. In essence, there seems to exist a magnetic field 
strength threshold above which the magnetic field acts to reduce the convection and below which 
the applied field actually enhances the convection. Our limited numerical experiments seem to 
suggest that in general a higher threshold value is needed for a higher Rayleigh number. 
Because of the limited time and resource allocated for this project, we were unable to quantify 
the effect in this report and also the possible effect of the other system parameters such as aspect 
ratios on the threshold value to obtain magnetic damping effect. This should be a very important 
issue that requires further study, considering the widespread use of magnetic damping in 
solidification processes. 

4. CONCLUDING REMARKS 

A fully 3-D numerical model was developed to represent magnetic damping of complex 
fluid flow, heat transfer and electromagnetic field distributions in a cylindrical melt cavity. The 
model was developed based on our in-house finite element code for the fluid flow, heat transfer 
and electromagnetic field calculations, which has been tested against benchmark test problems 
that are solved by other commercial codes as well as analytical solutions when available. The 
numerical model was tested against available numerical and experimental results reported for 
water. With the model so tested, various numerical simulations were carried out for the Sn- 
35.5% Pb melt. Numerical results showed that magnetic damping can be effectively applied to 
reduce turbulence and flow levels in the melt undergoing solidification and over a certain 
threshold value a higher magnetic field resulted in higher velocity reduction. It is found also that 
for a fully 3-D representation of the magnetic damping effects, the induced electric field in the 
melt does not vanish, as some researchers suggested, and must in general be included even for 
highly conducting melts such as molten metal and semiconductors. Also, for the study of the 
melt flow instability, a long enough simulation time has to be applied to ensure the final fluid 
flow circulation pattern. Moreover, our numerical results suggested that there seems to exist a 
threshold value of applied magnetic field above which magnetic damping becomes possible and 
below which the convection in the melt is actually enhanced. Because of the limited financial 
resource allocated for the project, were ere unable to carry out extensive study on this effect, 
which should warrant further theoretical and experimental study. In conclusion, the numerical 
model developed in this study may now be used to study the fundamentals governing the 
magnetic damping phenomena in the solidifying melts and should serve as a useful tool for 
exploring necessary design parameters magnetic damping experiments. 
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APPEDIX I: PROBLEM STATEMENT 


To simplify the calculations, magnetic field damping effects on the fluid flow in a 3-D 
cylindrical cavity as shown in Figure Al, instead of the whole solidification process as appears 
in Figure 1, are considered. 


B 



Figure Al. Schematic representation of the 3-D model for magnetic damping studies. 


The mathematical equations governing the fluid flow, temperature distribution and 
electromagnetic field distributions in the melt, along with appropriate boundary conditions, are 
written as follows. 

GOVERNING EQUATIONS 
Continuity Equation: 


V*u = 0 


( 1 ) 


Momentum Equation: 


N 


d\y 

— +u« Vu 

dt 


= -V/> + /tV 2 u - pg(5, (T - T 0 ) + a m (E + u x B) x B 


Energy Equation: 


f dT 

pc A + U.VT 

y p [ St 


\ 


= kW 2 T 


( 2 ) 


(3) 
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Electromagnetic Field Equation: 


J = cr m (E + u x B) 

(4) 

V*J = 0 

(5) 

Vx E = 0 

(6) 

E = -VO 

(7) 

V 2 0 = V • (B x u) = B • (V x u) 

(8) 


Thus for the problem under consideration, the Maxwell equations reduces to one equation 
for an unknown electric potential, that is, Eq. (8), and other relevant quantities can be calculated 
once the potential is known. This simplification comes from the fact that for metallic liquid, the 
induced magnetic field b is proportional to the magnetic Reynolds number and that the magnetic 
Reynolds number is much smaller than imposed magnetic field B 0 , thereby B=B 0 . This curl free 
field allows us to write E as the gradient of an electric potential O. This, combined with the fact 
that the current density is a divergence free field, results in V z O = B • (Vx u) . 

The boundary conditions for the problem represent the physical constraints needed to 
solve the problem. For the present problem, the no slip conditions are applied at the solid walls, 
side walls are thermally insulated, the top and bottom walls are at the fixed temperatures and the 
electric current density is zero on the solid walls as no current can escape from the system, 
because the walls are electrically insulated as well. Mathematically, these conditions are given 
as follows. 


Side Walls: 



u = 0 


<30 

dn 


Top Wall: 


<30 A 

T = T, u = 0 = 0 

£ dn 


t = t h 


dO 

77 


Bottom Wall: 


u = 0 


0 



To facilitate numerical computations, it is more convenient to non-dimensionalize the 
above equations and boundary conditions. To do that, we use the standard velocity, time, and 
length scales for natural convection in encloures [1] and define the following dimensionless 
variables. 


u*=u hla x’ = x/h t'=ta/h 2 

p=ph 2 !pa 2 T* =(T-T l )/(T h -T l ) 


B* =B/B 0 


0*=<D /aB 0 E'=Eh/aB 0 J'=Jh/a m aB 0 


With the above definitions, the governing equations and the boundary conditions may be re- 
written in non-dimensional forms. 

Continuity Equation: 

V* • u* = 0 (9) 


Momentum Equation: 

^-V + U* • V*u* =-V>* +/VV* 2 u* -RaPrT'e z +ffa 2 /V(-V* O* +u* xB*)xB* (10) 
d t 


Energy Equation: 


dT' 

~d7 


+ u* 


• v*r* 


* 2 * 

V T 


( 11 ) 


Electromagnetic Field Equation: 

V*V = V • (b* x u* ) = B* • (V* xu*) 


( 12 ) 


11 



Boundary conditions 


Side Walls: 

* 

II 

o 

u* =0 



dn 


dn 

Top Wall: 

T’ =0 

u' =0 





dn 

Bottom Wall: 

T' =1 

u* =0 



dn 


The non-dimensionalization process often gives rise to a set of dimensionless parameters 
that characterize the performance of the system under study. For our problems, these parameters 
are given in the following table. 


DIMENSIONLESS SYSTEM PARAMETERS 


Dimensionless Number 

Definition 

Rayleigh number 

Ra = ^ T »- T ‘ )h> 


va 

Prandtl number 

Pr = — 


a 

Hartman number 

Ha = £A — 


V P 

Thermal Diffusivity 

k 

a = 


p c p 

Kinematic Viscosity 

v = -tL 


P 
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APPENDIX II: FINITE ELEMENT FORMULATION 


1. Discretization and Interpolation: 



Figure A2. Finite element discretization of the 3-D cavity for magnetic damping of fluid flow studies. 

We started with discretizing the domain into a finite number of elements as shown in 
Figure A2. Over each element, the dependent variables u, T and O are interpolated by their 
nodal point values and the shape functions of <p(x,) and for P we use yAxi) as the interpolation 
functions. Then the velocity, pressure, temperature and electric potential fields are approximated 
by: 




(13) 

P(x,.,o = i/ r P(0 


(14) 

T(x i ,t) = fT(t) 


(15) 

o(x „o = # r <i>(0 

i,j =1, 2 , 3 

(16) 


where Uj, P, T and <t> are vectors of approximate values of element nodal points. 

Our calculations employ 8-node elements, as shown in Figure A3. The interpolation 
functions expressed in local coordinate are given by: 
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<l> = {N j ,N j ,---,N j ,--) T ;</ = 1 


(17) 


Nj =ia + ^Xl + Wy) + 0 + ffy) (18) 

o 

where ^ j ,r) j ,q j =± 1, j = 1,2,3,..., 8. Note that ^ is chosen one order lower than the 
interpolation function for the velocity field to satisfy the Babuska-Brezzi stability condition. 



Figure A3. 8-node element in local coordinate system. 


The interpolated velocity and other field variables can not satisfy the governing equations 
exactly. Substituting the interpolation equations into the governing equations, we get the 
residual functions Ri, R2, R3 and R4, which respectively represent error functions of the 
momentum, mass conservation, energy and electric potential equations. To simplify the 
calculations, the pressure may be treated using the penalty method, by which the penalty form of 
the continuity equation becomes, 


P = -—V • u (19) 

e 

where e is the penalty parameter. 

Substituting the interpolation functions into the governing eqations, one has the following 
residual functions for the field variables. 


f 1 (^,U„P,T,0) = R 1 

Momentum 

(20) 

f 2 (^,U,,P) = R 2 

Penalty Continuum 

(21) 
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f 3 ( < MJ l ,T) = R 3 

Energy 

(22) 

u,,o) = r 4 

Electric Potential 

(23) 


2. The Galerkin Weighted Residual Method: 

Finite element formulation of nonlinear equations, as discussed in the present case, is carried 
out using the Galerkin Weighted Residual method. By this method, the residuals of the 
differential equations are made orthogonal to independent weighting functions, and these 
weighting functions are chosen to be the same as the shape functions for the unknowns: 


[R 2 if/dV = 0 


(24) 

o 

II 

of 

i=l,3,4 

(25) 


Applying the Green theorem to the above equations and carrying out involved algebraic 
operations, one has the following set of matrix equations for the unknowns defined at the nodal 
points, 


MU+K(U)U + C P P + C B T = F(U) 

(26) 

M r T+K r (U)T = F r 

(27) 

-fM p P = C/U 

(28) 

K o O> + C 0 U = F 0 

(29) 


where now U is a global vector containing all the nodal values of u, v and w. P, T and O are 
global pressure, temperature and electric potential vectors respectively. M is the mass matrix, 
K(U)=N (TJl+lCto+K^ is the stiffness matrix, where is the viscous stiffness matrix, and N(U) 

is the advection stiffiiess matrix. C p is gradient matrix, K B the buoyancy matrix and F(U) the 
"force" vector which incorporates the Lorentz force term and the boundary conditions. In the 
temperature, penalty continuum and electric potential equations, the subscripts T, P and O refer 
to the corresponding sets of coefficients. Solving the pressure from continuity equation and 
substituting back into the momentum equation, one has the final form of the matrix 
representation of the governing equations in matrix form. 


M 

0 

o' 

u" 


'k(U)-c p m p -'c p t 

K B 

0 

u' 


F(U)~ 

0 

M r 

0 

f 

+ 

0 

k t 

0 

T 

= 

F t 

0 

0 

0 

i 


K 0 

0 

M 0 

O 


K . 


( 30 ) 
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where the coefficient matrices are defined 
M = M r = j <f><j> T dV ; 

k t =k* = fit.vfdv- 

K(U) = N(U) + K <m ,+K^; 

K,^ = -|(®xBx {^> T v^dV or K 
C P = [<f>{V^)dV- 
K Ji = [RaPrg^ T dV- 

F(U) = {F 1 (U),F 2 (U),F 3 (U)} r ; 

F r ; 


M p = | y/y/ T dV ; = J PrVcp • V</> T dV 

N(U) = N r (U)=j^u»V^F 
K r (U) = N r (U) + K r 

,ik = ~ £ fa ijk £ jlm B l B m $ U k dV 

C/ = | ¥ <yf)dv 

K* =-y>B*(Vx<fi T u)dV 

Fj (U) = - Jov> dS + | Ha 2 Pr(-V O xB) J dV 

F ® =-[^E<f >dS 


where is the same for all 3 velocity components and u k refers to the kth component of the 
velocity. 

The above matrix equations are solved using the successive substitution method and the 
time derivatives are approximated using the implicit finite difference scheme. Both variable time 
steps with automatic error tracking and time step automation and fixed time steps may be 
applied. The majority of the results presented in the report are calculated using the fixed time 
steps. 
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Table 1: Parameters used for calculations 


Parameter 

Value 

Diameter of cylinder (d) 

7xl0' 3 m 

Aspect ration (h/d) 

0.5, 1.0, 3.0 

Top wall temperature ( 7 / ) 

520 K 

Magnetic field Bo 

0.0~ 0.5 T 

Gravity (g ) 

9.8 m s' 2 
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Table 2 Parameters of liquid metal used for calculations 


Property ( Sn-35.5% Pb ) 

Value 

Density of the melt ( p ) 

1.068xl0 4 kg m' 3 

Thermal expansion coefficient (J3) 

1.23x1 O' 4 K' 1 

Thermal conductivity of the melt (k) 

14.9 Wra-'K' 1 

Specific Heat ( C p ) 

173.70 J Kg' K' 1 

Viscosity of the melt (ji) 

2.410xl0' 3 Nsm' 2 

Electronic conductivity of the melt (cr m ) 

3.402xl0 7 Qm'' 

Thermal diffusivity of the melt (a) 

8.032x1 0‘ 6 m 2 s' 1 

Kinematic viscosity of the melt ( v) 

2.257xl0‘ 7 m 2 s' 1 

Associated dimensionless numbers 

Value 

Prandtl number (Pr) 

0.0281 

Rayleigh number ( Ra ) 

4.0x10 3 ~ 8.0x10 5 

Hartmann number (Ha) 

0-1247.5 

Time and velocity scales used to 
normalize the results presented in all the 
figures 

Value 

Time scale (t) 

1.53 (sec) 

Velocity scale (uo) 

2.29 xlO' 3 (m/s) 
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Table 3 Properties of water used for calculations 


Property 

Value 

Thermal expansion coefficient (fi) 

0.227x1 O' 3 K" 1 

Thermal diffusivity of the melt (a) 

1.145xl0* 7 m 2 s' 1 

Kinematic viscosity of the melt ( v) 

0.963x1 O' 6 m 2 s' 1 

Associated dimensionless numbers 

Value 

Prandtl number (Pr) 

6.7 

Rayleigh number ( Ra ) 

2.88x10 3 ~ l.OxlO 5 
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Ra=2800 

h/d=0.5 



Figure 2. Comparison of our 3-D finite element model predictions for the velocity and 
temperature distributions (a) with the finite difference results (b) reported by Muller's research 
group (Muller et al., 1984) for natural convection of water (H 2 0) in a cylindrical container heated 
from below. The parameters used for calculations: h/d=0.5, Ra=2800 and Pr=6.7 (H 2 0). 
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Figure 3. Comparison of our 3-D finite element model predictions for the velocity and 
temperature distributions (a) with the finite difference results (b) reported by Muller's research 
group (Muller et al., 1984) for natural convection of water (H 2 0) in a cylindrical container heated 
from below. The parameters used for calculations: h/d=l, Ra=17500 and Pr=6.7 (H 2 0). 












Figure 4. 3-D view of temperature profile evolution for Sn-35.5% Pb melt, Pr=0.0281, 
Ra=4000 (T boltom = 660 K), h/d=0.5 at (a) t=0; (b) t=0.040; (c) t=0.4(X); and (d) t=2.000 
(=2(XX) time steps, steady state). The contour values are the same as shown in Figure 6. 




Figure 5. 3-D view of velocity field evolution for Sn-35.5% Pb melt, Pr=0.0281, 
Ra=40(X) (Tbooom = 660 K), h/d=0.5 at (a) t=0; (b) t=0.040; (c) t=0.400; and (d) 
t=2.0(X) (steady state). 





Figure 6. Transient development of flow fields and temperature profiles for Sn-35.5% Pb 
viewed in the x-z and y-z planes: Ra=4000 (Tbo«om = 660 K), Pr=0.0281, h/d=0.5; (a) and 
(b) t=0.20; (c) and (d) t=0.400; (e) and (f) t=1.000; (g) and (h) t=2.000 (steady state). 




Figure 7. 3-D view of temperature profile evolutation for Sn-35.5%Pb alloy 
Pr=0.0281, Ra=l() 4 (T telom =564 K), h/d=1.0: (a) t=0; (b) t=0.1(); (c) t=0.80; and 
t=2.20 (time steps=2200, steady state). The contour values are given in Figure 9. 







Figure 8. 3-D view of flow field evolutation for Sn-35.5%Pb Pr=().()281, Ra=10 4 , 
(Tbo«om=564 K), h/d=1.0: (a) t=0.02; (b) t=0.1; (c) t=0.8; and t=2.2 (steady state). 
















Figure 10. Fluid velocity and temperature profiles in the Sn-35.5% Pb melt, Pr=0.0281, 
Ra=10 ;> (Tbottom = 536 K), h/d=3.0, (time steps=1100, steady state): (a) 3-D view of fluid 
field; (b) and (c) x-z and y-z plane eut views of flow patten and temperature distribution. 
(Umax=0.53). 













Figure 13. Transient development of flow fields and temperature profiles in the Sn- 
35.5% Pb melt, Pr=(M)281, Ra=4(XX), h/d=0.5 with magnetic field B=0.05T applying at 
t=0.4; (a) and (h) t=0.4; (c) and (d) t=0.420; (e) and (f) t=0.5; and (g) and (h) t=0.8. 




Figure 14. Transient calculation for Pr=0.0281, Ra=10 4 , h/d=1.0 with magnetic Field 
B=0.()16T switched on at t=0.400: (a) and (b) 3-D view of fluid field and temperature 
profile; (c) and (d) the x-z and y-z plane cut views of flow pattern and temperature 
distribution (t=0.8, Umax=0.93). 




Figure 15. Transient calculation for Pr=0.0281, Ra=10 4 , h/d= 1.0 with magnetic field 
B=0.05T switched on at t=0.400: (a) and (b) 3-D view of fluid field and temperature 
profile; (c) and (d) the x-z and y-z plane cut views of flow pattern and temperature 
distribution ft=0.8. Umax=0.42). 






Figure 16. Velocity and temperature profile in the Sn-35.5% Pb melt, Pr=0.0281, 
Ra=10 5 , h/d=3.0 with magnetic field B=0.05T switched on at t=0.4: (a) 3-D view of 
fluid field; (b) and (c) x-z and y-z plane cut views of flow pattern and temperature 
distribution (t=0. 8, Umax=0.19). 
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Figure 17. Melt flow and temperature profile in the Sn-35.5% Pb liquid pool, 
Pr=0.0281, Ra=4(.XX), h/d=0.5 with magnetic field B=0.5T turned on at t=0.40: (a) 3- 
D view of fluid flow field profile; (b) and (c) the x-z and y-z plane cut views of flow 
pattern and temperature distribution (1=0.8, Umax=0.0l9). 





Figure 18. Fluid flow and temperature distribution in the Sn-35.5% Pb melt. 
Pr=O.0281, Ra=10 4 , h/d= 1 .0 with magnetic field B=0.5T switched on at 1=0.40: (a) 
and (b) 3-D view of fluid field and temperature profile; (c) and (d) x-z and y-z plane 
cut views of flow pattern and temperature distribution (t=0.8, Umax=0.019). 





Figure 19. Transient calculation for Pr=0.0281, Ra=10\ h/d=3.() with magnetic field 
B=0.5T switched on at t=0.40: (a) and (b) 3-D view of fluid field and temperature 
profile; (c) and (d) x-/. and y-/ plane cut views of flow pattern and temperature 
distribution (t=0.8. Umax=0.046). 







Figure 20. Dependency of the turbulent melt (Sn-35.5% Ph) flow on the applied 
magnetic field, Pr=0.0281, Ra=8xl() 5 (Tbottom =4027 K), h/d=I.O with different 
magnetic field applying at t=0.40: (a) and (b) B=0, Umax=l 72.50; (c) and (d) 
B=0.005T, Umax=96.59; (e) and (f) B=0.05T, Umax=38.79; (g) and (h) B=0.5T, 
Umax=1.57, all at t=0.8. 



Figure 21. Fluid flow and temperature distribution in the Sn-35.5% Pb melt in 3-D cavity 
without the presence of an applied magnetic field. According to Muller’s numerical simulations 
(run time (dimensionless) up to time steps=600), no symmetric flow patterns are possible. 
However, our calculations showed that if a longer time is used, a symmetric flow pattern starts to 
evolve and gradually attains a steady state. For the case shown, at time step=500, the computed 
results are similar to Muller’s; and a steady state condition seems to be reached at time 
step=2400. The results here seem to suggest that additional work might be needed to further 
clarify some of the conclusions made in earlier studies by Muller’s group (Muller 1988). 






Figure 22. Comparison of calculated results with and without the electric field term, 
Pr=0.0281, Ra=4(XX), h/d=0.5: (a) x-z plane view of flow field and temperature 
distribution, B=0, Umax=2.7790; (b) x-z plane view of velocity and electric potential 
field. B=0.005T, Umax=2.2803: (c) x-z plane view of flow field and temperature 
distribution with E set to 0 and B=0.005T, Umax=2.3531; (d) x-z plane view of velocity 
and electric potential field. B=0.05T, Umax=0.2254; (e) x-z plane view of flow field and 
temperature distribution with E set to 0 and B=0.05T, LJmax=0.3345. 









Figure 23. Dependency of the magnetic damping effects on the applied magnetic field strength. 
Note Utal with a weak magnetic field, the convection is actually enhanced, rite tlireshold value ol 
the applied magnetic field to achieve damping effects seems to be slightly larger than 0.01 6T. 



